#define RNSIZE 1000000

#include <iostream>
#include <fstream>
#include <math.h>
#include "probConvTableGenerator.h"

#define INDEXSIZE 14348907L
//define INDEXSIZE 348907L

using namespace std;



signed char mySign[INDEXSIZE];
signed char absoluteIndexToInequality[INDEXSIZE][15];



static inline double entropy(double a) {
  return -a*log2(a);
}

const double h_to_I[15][15]={ {0,0,0,0,0,0,0,0,0,0,0,0,-1,0,1},
						{0,0,0,0,0,0,0,0,0,0,0,0,0,-1,1},
						{0,0,0,0,0,0,0,0,0,0,0,-1,0,0,1},
						{0,0,0,0,0,0,0,0,0,0,-1,0,0,0,1},
						{0,0,0,0,0,0,-1,0,0,0,0,0,1,1,1},
						{0,0,0,0,0,0,0,0,0,-1,0,1,1,0,1},
						{0,0,0,0,0,-1,0,0,0,0,1,0,1,0,1},
						{0,0,0,0,0,0,0,0,-1,0,0,1,0,1,1},
						{0,0,0,0,0,0,0,-1,0,0,1,0,0,1,1},
						{0,0,0,0,-1,0,0,0,0,0,1,1,0,0,1},
						{0,-1,0,0,1,1,0,0,0,1,-1,-1,-1,0,1},
						{0,0,-1,0,0,1,1,1,0,0,-1,0,-1,-1,1},
						{0,0,0,-1,0,0,1,0,1,1,-0,-1,-1,-1,1},
						{-1,0,0,0,1,0,0,1,1,0,-1,-1,0,-1,1},
						{1,1,1,1,-1,-1,-1,-1,-1,-1,1,1,1,1,1}
						};

const int permutationTable[24][15]={	{1,2,4,3,5,7,6,9,8,10,11,13,12,14,15},
								{1,4,3,2,7,6,5,10,9,8,13,12,11,14,15},
								{4,2,3,1,9,10,7,8,5,6,11,12,14,13,15},
								{1,3,2,4,6,5,7,8,10,9,12,11,13,14,15},
								{3,2,1,4,8,6,10,5,9,7,11,14,13,12,15},
								{2,1,3,4,5,8,9,6,7,10,14,12,13,11,15},
								{1,3,4,2,6,7,5,10,8,9,12,13,11,14,15},
								{1,4,2,3,7,5,6,9,10,8,13,11,12,14,15},
								{3,2,4,1,8,10,6,9,5,7,11,13,14,12,15},
								{4,2,1,3,9,7,10,5,8,6,11,14,12,13,15},
								{2,4,3,1,9,8,5,10,7,6,13,12,14,11,15},
								{4,1,3,2,7,10,9,6,5,8,14,12,11,13,15},
								{3,1,2,4,6,8,10,5,7,9,14,11,13,12,15},
								{2,3,1,4,8,5,9,6,10,7,12,14,13,11,15},
								{2,1,4,3,5,9,8,7,6,10,14,13,12,11,15},
								{3,4,1,2,10,6,8,7,9,5,13,14,11,12,15},
								{4,3,2,1,10,9,7,8,6,5,12,11,14,13,15},
								{2,3,4,1,8,9,5,10,6,7,12,13,14,11,15},
								{3,1,4,2,6,10,8,7,5,9,14,13,11,12,15},
								{3,4,2,1,10,8,6,9,7,5,13,11,14,12,15},
								{2,4,1,3,9,5,8,7,10,6,13,14,12,11,15},
								{2,3,4,1,8,9,5,10,6,7,12,13,14,11,15},
								{4,1,2,3,7,9,10,5,6,8,14,11,12,13,15},
								{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15}
							};
				

int main(void)
{

	//generate the large indexedInquality 
	
	
	for(long i=0;i<INDEXSIZE;i++)
	{
		
		for(int j=0; j<=14;j++)
		{
			absoluteIndexToInequality[i][j]=(signed char)((i%((long)pow((double)3,(double)(j+1)))-i%((long)pow((double)3,(double)((j))))/(long)pow((double)3,(double)j))-1);
		}
		
		
		//mySign[i]=0;
		//~ if(i%((long)pow((double)10,(double)4))==0)
			//~ cout<<i<<endl;
	}
	cout<<"absoluteIndexToInequality initialized\n";
	
	//read in the random numbers from the file
	double rnBuffer[RNSIZE];

	ifstream rnIn( "random.data" );
	if( rnIn )
	{
		rnIn.read((char*)rnBuffer,RNSIZE*sizeof(double)/sizeof(char));
		rnIn.close();
		cout<<"random.data loaded in\n";
	}
	
	
	
	cout<<"first 16 random nubmers: \n";
	for(int i=0;i<16;i++)
	{
		cout<< rnBuffer[i] <<" ,";
	}
	cout<<endl;
	

	//only need to generate it once
	//probConvTableGenerator();
	
	indexedProb probConvTableIn;
	
	ifstream probConvTableInStream("probConvTable.data", ios::binary);
	if(probConvTableInStream)
	{
		probConvTableInStream.read((char*)&probConvTableIn, sizeof(indexedProb));
		probConvTableInStream.close();
		cout<<"probConvTable.data loaded in\n";
	}

	//cout<<"\n the table in\n";
	//printIndexedProb(probConvTableIn);
	long totalRandomSet=RNSIZE/16;
	
	double currentProbVector[16];
	double normalizeFactor=0;
	
	double hVector[15];
		//~ hVector[0]=h1;hVector[1]=h2;hVector[2]=h3;hVector[3]=h4;
		//~ hVector[4]=h12;hVector[5]=h23;hVector[6]=h34;hVector[7]=h13;hVector[8]=h14;hVector[9]=h24;
		//~ hVector[10]=h123;hVector[11]=h124;hVector[12]=h234;hVector[13]=h134;
		//~ hVector[14]=h1234;
	
	double &h1234=hVector[14];
	double &h123=hVector[10],&h124=hVector[11],&h134=hVector[13],&h234=hVector[12];
	double &h12=hVector[4],&h13=hVector[7],&h14=hVector[8],&h23=hVector[5],&h24=hVector[9],&h34=hVector[6];
	double &h1=hVector[0],&h2=hVector[1],&h3=hVector[2],&h4=hVector[3];
	double probability;
	
	double iVector[15];
	double permutedIvector[15];
	
	double inequalityValue=0;
	
	
	long ruledout=0, oldRuledOut=0;
	double minPositive;
	double maxNegative;
	minPositive=pow(10,-8);
	maxNegative=-pow(10,-8);
	
	ifstream prunedIn( "pruned.data" , ios::binary);
	if(prunedIn)
	{		
		prunedIn.read((char*)mySign,INDEXSIZE);
		prunedIn.close();
		cout<<"pruned.data loaded in\n";
	}
	
	//totalRandomSet
	for(long rnCounter=0; rnCounter<totalRandomSet; rnCounter++)
	{
		normalizeFactor=0;
		for(int i=0;i<16;i++)
		{
			normalizeFactor+=rnBuffer[rnCounter*16+i];
			currentProbVector[i]=rnBuffer[rnCounter*16+i];
		}
		//cout<<"normalizeFactor="<<normalizeFactor<<endl;
		
		for(int i=0;i<16;i++)
		{
			currentProbVector[i]/=normalizeFactor;
		}
		
		h1234=0;
		for(int i=0;i<16;i++)
		{
			h1234+=entropy(currentProbVector[i]);
		}
		//cout<<"h1234="<<h1234<<endl;
		//cout<<"hVector[14]="<<hVector[14]<<endl;
		
		h123=0;
		for(int i=0;i<8;i++)
		{
			h123+=entropy(currentProbVector[probConvTableIn.p123[i][0]]+currentProbVector[probConvTableIn.p123[i][1]]);
		}
		//cout<<"h123="<<h123<<endl;
		
		h124=0;
		for(int i=0;i<8;i++)
		{
			h124+=entropy(currentProbVector[probConvTableIn.p124[i][0]]+currentProbVector[probConvTableIn.p124[i][1]]);
		}
		//cout<<"h124="<<h124<<endl;
		
		h134=0;
		for(int i=0;i<8;i++)
		{
			h134+=entropy(currentProbVector[probConvTableIn.p134[i][0]]+currentProbVector[probConvTableIn.p134[i][1]]);
		}
		//cout<<"h134="<<h134<<endl;
		
		h234=0;
		for(int i=0;i<8;i++)
		{
			h234+=entropy(currentProbVector[probConvTableIn.p234[i][0]]+currentProbVector[probConvTableIn.p234[i][1]]);
		}
		//cout<<"h234="<<h234<<endl;
		
		h12=0;
		for(int i=0;i<4;i++)
		{
			probability=0;
			for(int j=0;j<4;j++){probability+=currentProbVector[probConvTableIn.p12[i][j]];}h12+=entropy(probability);
		}
		//cout<<"h12= "<<h12<<endl;
		
		h13=0;
		for(int i=0;i<4;i++)
		{
			probability=0;
			for(int j=0;j<4;j++){probability+=currentProbVector[probConvTableIn.p13[i][j]];}h13+=entropy(probability);
		}
		//cout<<"h13= "<<h13<<endl;

		
		h14=0;
		for(int i=0;i<4;i++)
		{
			probability=0;
			for(int j=0;j<4;j++){probability+=currentProbVector[probConvTableIn.p14[i][j]];}h14+=entropy(probability);
		}
		//cout<<"h14= "<<h14<<endl;
		
		h23=0;
		for(int i=0;i<4;i++)
		{
			probability=0;
			for(int j=0;j<4;j++){probability+=currentProbVector[probConvTableIn.p23[i][j]];}h23+=entropy(probability);
		}
		//cout<<"h23= "<<h23<<endl;
		
		h24=0;
		for(int i=0;i<4;i++)
		{
			probability=0;
			for(int j=0;j<4;j++){probability+=currentProbVector[probConvTableIn.p24[i][j]];}h24+=entropy(probability);
		}
		//cout<<"h24= "<<h24<<endl;
		
		
		h34=0;
		for(int i=0;i<4;i++)
		{
			probability=0;
			for(int j=0;j<4;j++){probability+=currentProbVector[probConvTableIn.p34[i][j]];}h34+=entropy(probability);
		}
		//cout<<"h34= "<<h34<<endl;
		
		h1=0;
		for(int i=0;i<2;i++)
		{
			probability=0; for(int j=0;j<8;j++){probability+=currentProbVector[probConvTableIn.p1[i][j]];}
			h1+=entropy(probability);
		}
		//cout<<"h1= "<<h1<<endl;
		
		h2=0;
		for(int i=0;i<2;i++)
		{
			probability=0; for(int j=0;j<8;j++){probability+=currentProbVector[probConvTableIn.p2[i][j]];}
			h2+=entropy(probability);
		}
		//cout<<"h2= "<<h2<<endl;
		
		h3=0;
		for(int i=0;i<2;i++)
		{
			probability=0; for(int j=0;j<8;j++){probability+=currentProbVector[probConvTableIn.p3[i][j]];}
			h3+=entropy(probability);
		}
		//cout<<"h3= "<<h3<<endl;
		
		h4=0;
		for(int i=0;i<2;i++)
		{
			probability=0; for(int j=0;j<8;j++){probability+=currentProbVector[probConvTableIn.p4[i][j]];}
			h4+=entropy(probability);
		}
		//cout<<"h4= "<<h4<<endl;
		
		for(int i=0;i<15;i++)
		{
			iVector[i]=0;
			for(int j=0;j<15;j++)
			{
				iVector[i]+=h_to_I[i][j]*hVector[j];
			}
		}
		
		// cout<<"the I vector: ";
		ruledout=0;
		for (int i=0;i<INDEXSIZE;i++)
		{
			
			inequalityValue=0;
			bool isRuledOut=false;
			for(int permutationIndex=0;permutationIndex<24;permutationIndex++)
			{
				for (int k=0;k<15;k++)
				{
					permutedIvector[k]=iVector[permutationTable[permutationIndex][k]];
				}
				
				for(int j=0;j<15;j++)
				{
					inequalityValue+=permutedIvector[j]*absoluteIndexToInequality[i][j];
				}
			
				//signed char a= mySign[i];
				if(mySign[i]==2)
				{
					isRuledOut=true;
					//cout<<"old ruledout="<<ruledout<<","<<"index="<<i<<endl;
				}
				else if(mySign[i]==0)
				{
					if(inequalityValue>minPositive){
						mySign[i]=1;
						}
					else if(inequalityValue<maxNegative){
						mySign[i]=-1;
						}
				}
				else
				{ 
					if((inequalityValue>minPositive)&&(mySign[i]==-1)){
						mySign[i]=2; isRuledOut=true;
						cout<<"new ruledout="<<ruledout<<","<<"index="<<i<<", at sample="<<rnCounter<<endl;
						}
					else if((inequalityValue<maxNegative)&&(mySign[i]==1)){
						mySign[i]=2;isRuledOut=true;
						cout<<"new ruledout="<<ruledout<<","<<"index="<<i<<", at sample="<<rnCounter<<endl;
						}
				}
			}
			if(isRuledOut){ruledout++;}
		}
	// cout<<endl;
	
		if(ruledout>oldRuledOut)
		{
			ofstream prunedOut( "pruned.data" , ios::binary);
			if(prunedOut)
			{		
				prunedOut.write((char*)mySign,INDEXSIZE);
				prunedOut.close();
				cout<<"pruned.data saved\n";
			}
		}
		cout<<"\n\n total ruledout="<<ruledout<<" ("<<((double)ruledout/(double)INDEXSIZE)<<" \%)"<<", "<<"delta= "<<ruledout-oldRuledOut<<" inequalities"<<", at sample "<<rnCounter<<endl<<endl;
		oldRuledOut=ruledout;
	}

	cout<<"Total ruledout="<<ruledout<<endl;
	return 0;
	
}


int probConvTableGenerator()
{
//buid entropy computation table
	
	indexedProb probConvTable;
	
	int tempIndex1;
	int tempIndex2;
	int currentAbsoluteIndex;
	
	for(int i1=0;i1<=1;i1++)
	{
		for(int i2=0;i2<=1;i2++)
		{
			for(int i3=0;i3<=1;i3++)
			{
				for(int i4=0;i4<=1;i4++)
				{
					currentAbsoluteIndex=i1*8+i2*4+i3*2+i4;
					
					//cout<<"currentAbsoluteIndex="<<currentAbsoluteIndex<<endl;
					
					//*******************************
					tempIndex1=i1*4+i2*2+i3;
					probConvTable.p123[tempIndex1][i4]=currentAbsoluteIndex;
					
					tempIndex1=i1*4+i2*2+i4;
					probConvTable.p124[tempIndex1][i3]=currentAbsoluteIndex;
					
					tempIndex1=i1*4+i3*2+i4;
					probConvTable.p134[tempIndex1][i2]=currentAbsoluteIndex;
					
					tempIndex1=i2*4+i3*2+i4;
					probConvTable.p234[tempIndex1][i1]=currentAbsoluteIndex;
					
					//*******************************
					tempIndex1=i1*2+i2*1;
					tempIndex2=i3*2+i4*1;
					probConvTable.p12[tempIndex1][tempIndex2]=currentAbsoluteIndex;
					
					tempIndex1=i1*2+i3*1;
					tempIndex2=i2*2+i4*1;
					probConvTable.p13[tempIndex1][tempIndex2]=currentAbsoluteIndex;
					
					tempIndex1=i1*2+i4*1;
					tempIndex2=i2*2+i3*1;
					probConvTable.p14[tempIndex1][tempIndex2]=currentAbsoluteIndex;
					
					tempIndex1=i2*2+i3*1;
					tempIndex2=i1*2+i4*1;
					probConvTable.p23[tempIndex1][tempIndex2]=currentAbsoluteIndex;
					
					tempIndex1=i2*2+i4*1;
					tempIndex2=i1*2+i3*1;
					probConvTable.p24[tempIndex1][tempIndex2]=currentAbsoluteIndex;
					
					tempIndex1=i3*2+i4*1;
					tempIndex2=i1*2+i2*1;
					probConvTable.p34[tempIndex1][tempIndex2]=currentAbsoluteIndex;
					
					
					//*******************************
					tempIndex2=i2*4+i3*2+i4;
					probConvTable.p1[i1][tempIndex2]=currentAbsoluteIndex;
					
					tempIndex2=i1*4+i3*2+i4;
					probConvTable.p2[i2][tempIndex2]=currentAbsoluteIndex;
					
					tempIndex2=i1*4+i2*2+i4;
					probConvTable.p3[i3][tempIndex2]=currentAbsoluteIndex;
					
					tempIndex2=i1*4+i2*2+i3;
					probConvTable.p4[i4][tempIndex2]=currentAbsoluteIndex;
					
				}
				
			}
			
		}
	}
	
	cout<<"\n the table out\n";
	//printIndexedProb(probConvTable);
	
	ofstream probConvTableOutStream( "probConvTable.data", ios::binary );
	if( probConvTableOutStream )
	{
		probConvTableOutStream.write((char*)&probConvTable, sizeof(indexedProb));
		probConvTableOutStream.close();
	}
	
	indexedProb probConvTableIn;
	
	ifstream probConvTableInStream("probConvTable.data", ios::binary);
	if(probConvTableInStream)
	{
		probConvTableInStream.read((char*)&probConvTableIn, sizeof(indexedProb));
		probConvTableInStream.close();
	}
	
	cout<<"\n the table in\n";
	//printIndexedProb(probConvTableIn);
	
	return 0;
}


int printIndexedProb(struct indexedProb t)
{
	cout<<"\n p123: ";
	for(int i=0;i<8;i++)
	{
		for(int j=0;j<2;j++)
		{
			cout<<t.p123[i][j]<<",";
		}
		cout<<" ; ";
	}
	cout<<"\n";
	
	cout<<"\n p124: ";
	for(int i=0;i<8;i++)
	{
		for(int j=0;j<2;j++)
		{
			cout<<t.p124[i][j]<<",";
		}
		cout<<" ; ";
	}
	cout<<"\n";
	
	
	cout<<"\n p234: ";
	for(int i=0;i<8;i++)
	{
		for(int j=0;j<2;j++)
		{
			cout<<t.p234[i][j]<<",";
		}
		cout<<" ; ";
	}
	cout<<"\n";
	
	cout<<"\n p134: ";
	for(int i=0;i<8;i++)
	{
		for(int j=0;j<2;j++)
		{
			cout<<t.p134[i][j]<<",";
		}
		cout<<" ; ";
	}
	cout<<"\n";
	
	cout<<"\n p12: ";
	for(int i=0;i<4;i++)
	{
		for(int j=0;j<4;j++)
		{
			cout<<t.p12[i][j]<<",";
		}
		cout<<" ; ";
	}
	cout<<"\n";
	
	cout<<"\n p13: ";
	for(int i=0;i<4;i++)
	{
		for(int j=0;j<4;j++)
		{
			cout<<t.p13[i][j]<<",";
		}
		cout<<" ; ";
	}
	cout<<"\n";
	
	cout<<"\n p14: ";
	for(int i=0;i<4;i++)
	{
		for(int j=0;j<4;j++)
		{
			cout<<t.p14[i][j]<<",";
		}
		cout<<" ; ";
	}
	cout<<"\n";
	
	cout<<"\n p23: ";
	for(int i=0;i<4;i++)
	{
		for(int j=0;j<4;j++)
		{
			cout<<t.p23[i][j]<<",";
		}
		cout<<" ; ";
	}
	cout<<"\n";
	
	cout<<"\n p24: ";
	for(int i=0;i<4;i++)
	{
		for(int j=0;j<4;j++)
		{
			cout<<t.p24[i][j]<<",";
		}
		cout<<" ; ";
	}
	cout<<"\n";
	
	cout<<"\n p34: ";
	for(int i=0;i<4;i++)
	{
		for(int j=0;j<4;j++)
		{
			cout<<t.p34[i][j]<<",";
		}
		cout<<" ; ";
	}
	cout<<"\n";
	
	cout<<"\n p1: ";
	for(int i=0;i<2;i++)
	{
		for(int j=0;j<8;j++)
		{
			cout<<t.p1[i][j]<<",";
		}
		cout<<" ; ";
	}
	cout<<"\n";
	
	
	cout<<"\n p2: ";
	for(int i=0;i<2;i++)
	{
		for(int j=0;j<8;j++)
		{
			cout<<t.p2[i][j]<<",";
		}
		cout<<" ; ";
	}
	cout<<"\n";
	
	
	cout<<"\n p3: ";
	for(int i=0;i<2;i++)
	{
		for(int j=0;j<8;j++)
		{
			cout<<t.p3[i][j]<<",";
		}
		cout<<" ; ";
	}
	cout<<"\n";
	
	
	cout<<"\n p4: ";
	for(int i=0;i<2;i++)
	{
		for(int j=0;j<8;j++)
		{
			cout<<t.p4[i][j]<<",";
		}
		cout<<" ; ";
	}
	cout<<"\n";

	return 0;
}
